Physics of Solid and Liquid Alkali Halide Surfaces Near the Melting Point 



in 

O 

o 

O 

Q 



-a 

G 

O 

o 



> 

(N 

oo 

(N 

in 
o 

-I— > 
c3 



c 

o 
o 



X 



T. Zykova-Timan, 1 D. Ceresoli, 1 U. Tartaglino, 1 ' 2 and E. Tosatti 1 ' 3 

1 International School for Advanced Studies (SISSA), 
and INFM DEMOCRITOS National Simulation Center, 
via Beirut 2-4, 1-34014 Trieste, Italy 
2 IFF, FZ-Jiilich, 52425 Julich, Germany 
'International Centre for Theoretical Physics (ICTP), P.O.Box 586, 1-34014 Trieste, Italy 

This paper presents a broad theoretical and simulation study of the high temperature behavior of 
crystalline alkali halide surfaces typified by NaCl(lOO), of the liquid NaCl surface near freezing, and 
of the very unusual partial wetting of the solid surface by the melt. Simulations are conducted using 
two-body rigid ion BMHFT potentials, with full treatment of long-range Coulomb forces. After a 
preliminary check of the description of bulk NaCl provided by these potentials, which seems generally 
good even at the melting point, we carry out a new investigation of solid and liquid surfaces. Solid 
NaCl(lOO) is found in this model to be very anharmonic and yet exceptionally stable when hot. It 
is predicted by a thermodynamic integration calculation of the surface free energy that NaCl(lOO) 
should be a well ordered, non-melting surface, metastable even well above the melting point. By 
contrast, the simulated liquid NaCl surface is found to exhibit large thermal fluctuations and no 
layering order. In spite of that, it is shown to possess a relatively large surface free energy. The 
latter is traced to a surface entropy deficit, reflecting some kind of surface short range order. Finally, 
the solid-liquid interface free energy is derived through Young's equation from direct simulation of 
partial wetting of NaCl(lOO) by a liquid droplet. It is concluded that three elements, namely the 
exceptional anharmonic stability of the solid (100) surface, the molecular short range order at the 
liquid surface, and the costly solid liquid interface, all conspire to cause the anomalously poor 
wetting of the (100) surface by its own melt in the BMHFT model of NaCl - and most likely also 
in real alkali halide surfaces. 

PACS numbers: 68.35.Rh, 68.35.Md, 68.45.Gd, 82.65.Dp, 68.10.Cr, 61.50.Jr 



I. INTRODUCTION 

Attention is increasing toward adhesion and wetting, 
the structure and physics of solid-liquid interfaces espe- 
cially at high temperatures, and the structure of liquid 
surfaces, particularly of complex and molecular systems. 
In order to gain more insight into these problems, there 
is a strong need for good case studies, to use as well- 
understood starting points. 

One easy starting point is to study the contact of a 
liquid with its own solid, a clear situation where there 
will be no ambiguity of physical description, no uncer- 
tainty in chemical composition, no segregation phenom- 
ena, all of them complications present in the study of 
contact between different substances. Contact of a liq- 
uid with the surface of its own solid usually material- 
izes spontaneously at high temperature. Most solid sur- 
faces are known to wet themselves spontaneously with an 
atomically thin film of melt, when their temperature T 
is brought close enough to the melting point T m of the 
solid. The phenomenon whereby the thickness l(T) of 
the liquid film diverges continuously (and critically) as 
T — > T m , is commonly referred to as (complete) surface 
melting [111. 

Surface melting is indeed a very natural thing to hap- 
pen, because it corresponds to complete wetting of a solid 
substrate by the same identical substance, only in liquid 
form. Another name for surface melting, better suited 
for the fluid community 0, , could be complete inter- 
facial wetting. Due to surface melting, a solid with free 



surfaces cannot generally be overheated above its ther- 
modynamical bulk T m . The free energy barrier for the 
passage from solid to liquid generally requires nucleation 
and implies hysteresis. Although this is not generally 
observable for reasons exposed below, it should be theo- 
retically possible to overheat a solid, at least insofar as 
one can exclude the presence of defects that could act as 
nucleation centers. 

In the melting phase transition, the solid surface itself 
represents a (nearly) ubiquitous defect that continuously 
nucleates the liquid, thus usually preventing overheating 
of the solid. A generic crystal surface simply does not re- 
main solid at high temperature, and spontaneously wets 
itself with a liquid film of increasing thickness very close 
to the melting point. The thermodynamical condition for 
surface melting, or complete interfacial wetting, to occur, 
is that the solid-vapor (SV) interface should turn itself 
spontaneously into the solid-liquid (SL) plus liquid- vapor 
(LV) interface pair, or: 



7s v = 7sl + TlVi 



(1) 



where the 7's are interface free energies at the triple 
point. 

There are a number of known exceptions to this be- 
havior. On some solid surfaces the liquid film does be- 
gin to form upon heating, but its thickness levels off to 
a finite value instead of diverging as T — > T m (so-called 
blocked or incomplete |5| surface melting). More remark- 
ably, some other solid surfaces remain dry and fully crys- 
talline up to the bulk triple point. This surface non- 
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melting phenomenon, originally discovered in molecular 
dynamics simulations of Au(lll) [al and independently 
observed experimentally in Pb(lll) 0, takes place at the 
close-packed faces of several metals, such as Al(lll) 

Thcrmodynamically, surface non-melting will occur if 
there is a free energy loss in converting the SV interface 
into the SL plus LV interfaces pair, namely 

7sv < 7sl + 7lv- (2) 

In that case, the liquid will wet the solid at best in- 
completely. A liquid droplet deposited in full equilibrium 
on the solid surface will not spontaneously spread. It will 
settle instead in a metastable partial wetting geometry 
such as that of Fig. ^ The metastable droplet is char- 
acterized by two angles 9sl,0lv, and by two curvature 
radii (not shown) Rsl, Rlv of the solid-liquid and liquid- 
vapor interfaces, approximately obeying the generalized 
Young equations 0: 

7sv = 7s l cos 9 S l + 7lv cos 9 LV (3) 
R L vsm9 LV = RsLsin9sL (4) 

The temporary settling of a metastable liquid droplet on 
the surface of the same solid substance, schematically de- 
picted in Fig.^and discussed by Nozieres 0, was demon- 
strated in simulation in Ref. |lCI for A1/A1(111), but not 
verified experimentally yet [ll|. 

Here we wish to move on from elemental to molec- 
ular systems. Alkali halides are our natural prototype 
choice, because they represent a well defined class of sub- 
stances whose liquids do not wet their own solid, and 
because they have otherwise long been studied experi- 
mentally and theoretically. Molten salts and their sur- 
faces were extensively investigated by macroscopic tech- 
niques (3 ■ A partial wetting angle 9lv— 48° is known 
for liquid NaCl on NaCl(lOO) at the melting point T m 
= 1074 K (and similar results also hold for other al- 
kali halides) 0, [w|. Such a large angle underlines a 
strikingly poor wetting of the liquid onto its own solid, 
large by comparison with other known cases. For ex- 
ample, 9lv ^15° is observed for liquid Pb/Pb(lll) [la, 
or ~ 19° is obtained by simulation of Au/Au(lll) [lOj . 
Should liquid surface layering be, as was the case in 
metals, the culprit in NaCl too, the layering magnitude 
and its effect should indeed be exceptionally strong, and 
amenable to experimental verification. Several simula- 
tions on alkali halide li quid s suggest instead that there is 
no layering whatsoever [l(| UJ$ • A second possible mech- 
anism leading to surface non-melting may arise from van 
der Waals forces. Whenever the melt is optically more 
dense than the solid, the so-called Hamaker constant H 
governing the effective SL-LV interface interaction H/l 2 
may turn negative. The resulting attraction impedes 
complete melting, as is the case for valence semiconduc- 
tors such as silicon 0. However, liquid NaCl is 26% less 
dense than the solid, and here the Hamaker constant is 
certainly positive |l8j . 

That leaves the question of explaining the exception- 
ally poor wetting of NaCl(lOO) by liquid NaCl completely 



open. In order to shed light on this question and on 
the underlying physics, we undertook extensive simula- 
tion studies of the molten NaCl surface, and also of the 
NaCl(lOO) solid surface, at and close to the melting point. 

The plan of the rest of this paper is as follows. We 
will introduce first the simulation methods and the po- 
tentials used. Subsequently we will proceed to a careful 
characterization of all bulk properties of this model. The 
bulk solid and the bulk liquid will be simulated, and the 
results compared to experimental data and to previous 
theory and simulations. The bulk zero-pressure melting 
temperature, very close to the triple point temperature, 
will be extracted for this model potential by direct sim- 
ulation of the solid-liquid coexistence. Next, extensive 
slab simulations will be used to obtain a quantitative de- 
scription of the solid (100) surface near and above bulk 
melting, and of the liquid surface in the same temper- 
ature range. The solid surface free energy jsv will be 
calculated by thermodynamic integration. That of the 
liquid surface jlv will be calculated by means of the 
Kirkwood-Buff virial formula. The two will be compared 
and discussed. An effective harmonic calculation will be 
implemented showing the importance of anharmonicity 
in stabilizing the high temperature solid surface. A mod- 
ified calculation of the liquid surface free energy will be 
introduced to understand the poor temperature depen- 
dence of JlVi eventually explained in terms of surface 
short range order. Finally the solid-liquid interface free 
energy jsl will be calculated via Young's equation from 
a direct simulation of partial wetting of a liquid droplet 
on the solid surface at the melting point. The result- 
ing large value of 7sl is connected to the large density 
jump. In the concluding discussion, it is argued that all 
three separate physical mechanisms that conspire to give 
rise to poor wetting, by lowering jsv and simultaneously 
raising 7/,y and jsl eventually stem from charge order 
and charge neutrality of this ionic system. 

II. HAMILTONIAN, AND SIMULATION 
METHOD 

NaCl was modeled by the classic pairwise Born-Mayer- 
Huggins-Fumi-Tosi (BMHFT) rigid ion potential 19]: 

v ( r ij) = + A ap ex.p(-Brij) ^ g^. (5) 

Here a and j3 stand for either + or — , Z a and Zp are the 
ionic charges (+1 for Na and —1 for CI), the first term is 
the Coulomb interaction energy, the second is the short- 
range Pauli repulsion, and last two terms are induced 
dipole-dipole and dipole-quadrupole van der Waals in- 
teractions. The values of the parameters are reported in 
Table |H| 

Polarization forces 01 > though not negligible for the 
quantitative description of molecules, low-density ionic 
fluids, and other properties, do not appear to be crucial 
in the present context. They are neglected in order to 



FIG. 1: Sketch of the liquid drop partially wetting its own solid, showing the balance of the forces acting at the interfaces. 





Na-Na 


Cl-Cl 


Na-Cl 


A (eV) 


424.097 


3488.998 


1256.31 


B (A- 1 ) 


3.1545 


3.1545 


3.1545 


C (eV A 6 ) 


1.05 


72.5 


7.0 


D (eV A 8 ) 


0.499 


145.427 


8.676 



TABLE I: Parameters of Born-Mayer-Huggins-Fumi-Tosi po- 
tential for NaCl. 



concentrate on the BMHFT model, whose greater sim- 
plicity allows a much more extensive computational ex- 
ploration of the wetting properties [2(j • 

Bulk systems were studied at constant volume with 
cubic simulation cells comprising up to 10000, but more 
typically 3000-^5000, NaCl molecular units. Surfaces 
were studied with periodically repeated slabs - consisting 
of 12 -r 24 planes - separated by about 80 A of vacuum. 
The long range Coulomb interaction between ions was 
treated in full. We implemented a 3D Ewald summa- 
tion with repeated slabs and periodic boundary condi- 
tions along the direction z orthogonal to the surfaces as 
well as parallel to the (cc, y) surface plane. The alterna- 
tive choice of a single slab with 2D Ewald summation, 
physically more natural, was discarded as computation- 
ally more demanding [2l| . In our repeated slab 3D Ewald 
scheme, the large vacuum thickness is required in or- 
der to reduce spurious electrostatic couplings of instanta- 
neous fluctuating dipoles in liquid surfaces facing one an- 
other across the vacuum gap. Moreover, the Ewald sums 
were carried out for better convergence with conduct- 
ing ("tin foil") boundary conditions in the in-plane di- 
rections and with insulating or vacuum boundary condi- 
tions in the direction normal to the surface (see Ref. plf'). 
We performed both microcanonical and canonical simu- 
lations. In canonical runs, temperature was controlled 
by velocity-rescaling. Despite the size and time limita- 
tions imposed by long range forces, great care was taken 
to run simulations long enough for a clear equilibration, 
typically 100-^300 ps at T m , but longer when required. 
Finally, in all simulations involving the solid, whether as 
a bulk crystal, or at solid-liquid coexistence, or as a solid 
slab, we adjusted the cell size at each temperature so as 
to enforce the theoretical equilibrium lattice parameter 
independently computed as the one that gave a vanishing 
bulk stress (see Sec. IIII A"|) . 



III. BULK PROPERTIES 

A. Bulk Solid NaCl 

First of all, we verified how faithfully the bulk prop- 
erties of solid and liquid NaCl are reproduced by the 
BMHFT potentials. Thermal expansion of the solid is 
dictated by the increase of equilibrium lattice spacing 
at zero pressure as a function of temperature. We per- 
formed several simulations at constant temperature and 
volume, and computed the equilibrium lattice spacing 
by seeking the cell size that yielded vanishing pressure at 
each given temperature. The overall interpolated result 
is shown in Fig. [21 where it is compared with experi- 
mental data from various sources |22j . The temperature 
dependent cell size obtained in this manner was subse- 
quently enforced in all subsequent simulations involving 
solid NaCl, as mentioned earlier. The 298 K calculated 
equilibrium lattice spacing and linear expansion coeffi- 
cient a = {W^dV/dT were 5.683 A and 40.5-10~ 6 K" 1 
respectively (5.635 A and 38.3-10~ 6 K^ 1 are the experi- 
mental values [22|L 
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FIG. 2: NaCl volume expansion vs temperature. Both the 
theoretical (filled circles) melting temperature and change of 
the volume at T m are very similar to the experimental ones 
(empty triangles) . 

At higher temperature, notably above 600 K, the an- 
harmonicity is seen to get stronger, and expansion be- 
comes somewhat nonlinear. It should be noted that our 
procedure continues to work, and is in fact very accurate 
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at these higher temperatures. The simulated bulk solid 
remains locally stable and does not spontaneously melt 
(at least for our cell sizes and for times within 200 ps) 
until a maximum temperature T s ~ 1305 K. This maxi- 
mum bulk metastability temperature, necessarily higher 
than the ordinary melting temperature T m , approximates 
the "spinodal" temperature of bulk NaCl, defined as the 
point where the solid phase ceases to be a local free en- 
ergy minimum (for example by losing its mechanical sta- 
bility). As we will see further below, this bulk spinodal 
temperature is indeed predicted well above the melting 
temperature, precisely T s ~ T m + 240 K in the BMHFT 
model. It would be interesting to check this predisction 
by studying, e.g., spontaneous nucleation of the liquid in- 
side the solid, locally heated by for example two crossed 
laser beams. 

Root mean square displacements (Ar 2 ) 1 / 2 (RMSD) of 
Na + and Cl _ ions were extracted from the high tempera- 
ture bulk NaCl solid simulations. As shown in Fig. the 
value at 1066 K is 0.61 A for Na and 0.59 A for CI, com- 
parable with experimental estimates of 0.5 A and 0.48 A 
at the melting point |2^]. The corresponding calculated 
Lindemann ratios at 1066 K 6 — (Ar 2 ) 1 / 2 /a, (a is the in- 
teratomic distance) are 20% and 22% , compared with ex- 
perimental estimates at T m obtained from Debye- Waller 
factors of 17% and 20% |2^|. Simulation appears to 
slightly overestimate the thermal vibration amplitudes, 
but we note that the uncertainties in the experimen- 
tal procedure where RMSDs were extracted seem much 
larger than this discrepancy. In any case the RMSD val- 
ues are very much larger than the typical values between 
10% and 14% of the Lindemann ratio for most solids at 
T m . This large overshoot of the alkali halide bulk Linde- 
mann ratios can be rationalized by noting that whereas 
large thermal vibrations may very effectively destabilize 
atoms inside its bulk solid cage when interatomic forces 
are short ranged, they will much less effectively do so 
when forces are long range as is the case in strongly ionic 
solids. Simulations show that high temperature NaCl is 
in the BMHFT model a strongly vibrating, strongly an- 
harmonic, and yet unusually stable solid. As we shall 
show later, the same is true of the NaCl(100) surface. 

The vibrational spectral properties of warm NaCl can 
also be easily extracted, for later use, by Fourier trans- 
forming the velocity-velocity correlation functions: 



1 r°° 

^Na(w) = — / dt B Wt <V Na (t) • VNa(0)} , 
#K JO 



(6) 



and similarly for CI. The two spectral densities V^ a (ui) 
and Vci(w) extracted at 300 K and at 1066 K are shown 
m the Fig. H 

Generally, we note that our simulated vibrational spec- 
tra do not compare well in the details - even if not unrea- 
sonably in their gross features - with the more realistic 
spectral densities of Ref. |24| taken from the literature, 
and corresponding to very accurate shell model fits to ex- 
perimental phonon spectra (Fig.[3J). The rigid ion model 
is notoriously too crude to reproduce fine features as the 
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FIG. 3: The instantaneous displacements of the Na + (a) and 
CI - (b) ions in simulated bulk NaCl at 1066 K (bulk melting 
point of the Tosi-Fumi model). The displacement distribu- 
tions are well fit by gaussian distributions whose widths are 
indicated. The much smaller widths predicted by the Linde- 
mann melting criterion (Ar 2 ) 1//2 ~ 0.14 x a, where a is the 
interatomic distance, are also shown. 



detailed vibrational spectra, heavily affected by non rigid 
ion effects [2^ ■ Nonetheless it is not unreasonable to be- 
lieve that the overall temperature evolution of spectral 
densities can still be taken as representative of the real 
situation. The high temperature spectra show a consid- 
erable anharmonic softening and broadening relative to 
the room temperature ones. 
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FIG. 4: Vibrational density of states of bulk solid NaCl for (a) 
Na + and (b) Cl~ ions. Comparison between 300 K (solid line) 
and 1000 K (dashed line) indicates a considerable vibrational 
softening. 
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FIG. 5: Vibrational density of states of solid NaCl at low 
temperatures projected on (a) Na + and (b) CP ions. From 
Ref. 24] 
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B. Bulk Liquid NaCl 

Simulations of bulk liquid NaCl were conducted in 
close analogy to those of solid NaCl. In particular, the 
cell volume was adjusted isotropically to ensure vanishing 
pressure at every temperature. Our overall description of 
the liquid is very similar to that provided by many long- 
standing studies [IE 13 HE HI ■ The calculated liquid 
molecular volume is shown along with experimental val- 
ues in Fig. [21 The calculated volume expansion of 27% at 
melting compares very well with 26% from experiment. 
The internal structure of the liquid, notably the Na-Cl, 
Na-Na and Cl-Cl radial distribution functions: 

g(r)=p- 1 (Y / 6(r-R ij )y (7) 

quantities that are long well known, are well reproduced 
(Fig. HJ) by the simulation. They exhibit the charge cor- 
relation typical of molten salts, the first peak of <? H (r) 

giving a typical Na-Cl distance in the liquid of 2.6 A, 
about 10% shorter than in the solid. 




FIG. 6: The partial radial distribution functions at T m = 1066 
K. 



From the pair correlation functions, we may calculate 
the bulk liquid coordination number: 

N±=pf 4irr 2 g + ^(r)dr, (8) 
Jo 

where the upper integration limit corresponds to the first 
local (?-| (r) minimum, r m = 4.0 A. We find a coordina- 
tion number N = 4.6 in liquid NaCl at 1066 K, to be 
compared with the much higher N = 6 coordination in 
the defect-free bulk solid. Experimental estimates of the 
liquid coordination number vary in a broad range from 
N = 4.7 [28j to N = 5.8 |2flj. According to our simulation 
an ion in the liquid is surrounded by a cage consisting of 
only 4.6 ions of the opposite charge, albeit considerably 



closer than in the solid. In order to better understand 
the nature of the effective ion cage in the liquid we stud- 
ied the Na-Cl-Na and Cl-Na-Cl angular distribution in 
the bulk liquid: 

where |Ry|,|Rj fe | < r m . The result shown in Fig. [7| 
exhibits a main peak at 90° and only a weaker shoul- 
der around 150°). Similar results had been obtained by 
Amini 30] . The dominant 90° peak indicates in partic- 
ular that the local ionic cage surrounding an ion of the 
opposite charge is, even in the liquid, still roughly cubic, 
or more precisely octahedral, as in the solid. The octa- 
hedral cage linear size is smaller by about 10% than that 
of the solid. The 4.6 surrounding ions are distributed on 
the 6 corners of this cage, which therefore contains on 
average about 1.4 vacancies. 




9,(degree) 

FIG. 7: The angular distribution function of simulated liquid 
NaCl at T m . 



The diffusion coefficient: 




was also computed during simulations in order to check 
the general quality of description of the liquid. We found 
D = 10.53 • 10~ 5 cm 2 /s at T = 1300 K, which compares 
well with experimental values like D = 8.6 • 10 -5 cm 2 /s 
at T = 1121 K [3]]], as well as with those from previ- 
ous MD simulations, such as D = 9.5 • 10 ~ 5 cm 2 /s at 
T= 1267 K[27|. 

Based on all the above we conclude that the BMHFT 
potential description of bulk liquid NaCl is on the whole 
very good. 

C. Bulk Melting Temperature 

The melting temperature of the BMHFT model of bulk 
NaCl is obtained directly from simulations of solid-liquid 



6 



coexistence at zero pressure 32]. At constant tempera- 
ture, the interface in a bulk system which is roughly half 
solid and half melted (actually there are two SL inter- 
faces because of periodic boundary conditions) will drift 
with simulation time one way or the other so long as 
the system is away from the melting temperature T m , 
and will only remain stationary at T — T m . At constant 
energy, the average temperature (T) will instead slowly 
drift toward the melting temperature T m |3l, . 

We started with a crystalline bulk made up of 2880 
molecular units, a geometry comprising 6x6x20 conven- 
tional 4-molecule cubic NaCl cells along (x, y, z) respec- 
tively. Periodic boundary conditions (PBC) were as- 
sumed in all directions, and the volume was constant 
during each simulation. However, the cell size was ad- 
justed with temperature, so as to enforce zero stress at 
each temperature, as detailed below. After equilibration 
in proximity of the presumed melting point (T rs 1100 K) 
about one half of the 20 layers were melted by bringing 
them selectively at a higher temperature, while the re- 
maining atoms in the solid phase are kept fixed on their 
positions. The liquid and the solid halves initially out of 
equilibrium were subsequently let evolve to their mutual 
equilibrium. The in-plane cell size, and with it the in- 
plane solid lattice spacing were kept fixed at their solid 
equilibrium value previously established for that temper- 
ature as shown in Fig. [21 of Subsec. rill Al The linear cell 
size in the z-direction perpendicular to the solid-liquid 
interfaces was subsequently adjusted after each equilibra- 
tion so as to cancel the uniaxial stress normal to the inter- 
faces, compensating in particular the solid-liquid volume 
expansion. When the system was finally equilibrated by 
means of canonical MD runs (see Fig. |SJ, we observed 
the anticipated drift of the solid-liquid interfaces in op- 
posite directions at 1050 K and 1070 K, indicating that 
the melting temperature must fall in between (Fig. Et). 

To sharpen up our estimate of T m we then carried 
out delicately selected microcanonical simulation runs 
between E = —7.2 and —7.25 eV/molecule. Our final 
result was T m = 1066 ± 20 K (Fig. EJd). The conserva- 
tively large error bar quoted is mainly due to a residual 
pressure uncertainty of ~0.5 kbar. The real uncertainty 
is probably smaller, since we also found later that systems 
with free surfaces, and thus with a pressure more accu- 
rately close to zero, melt in bulk within 5 K of 1066 K. 
We also estimate that the error of T m due to other fac- 
tors, including small system size, and fluctuations, to be 
smaller than 20 K. For instance, the melting tempera- 
ture of an unrealistically small system made of 6x6x4 
molecules was found to be ~ 1050 K. Finally our value 
of the melting temperature, obtained in a totally unbi- 
ased manner |35j , is in essentially perfect agreement with 
T m = 1064 ± 14 K independently obtained by thermo- 
dynamic integration by Frenkel's group |36l ] . 




FIG. 8: Simulated coexistence of bulk solid and liquid. Peri- 
odic boundary conditions on all sides. 
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Simulation 


Experiment 


T m (K) 


1066±20 


1074 


AV 


27% 


26% 


L (eV/molecule) 


0.29 


0.29 


AS m (fc fl ) 


6.32 


6.38 


dP/dT (kbar/K) 


0.0311 


0.0357 [37] 


a (10" 6 K" 1 ) 


40.5 


38.3 


RMSD (A) 


0.60 


0.49 


S 


20-24% 


17-20% [23] 



TABLE II: High temperature properties of NaCl. T m is the 
melting temperature; AV is the volume jump at the melt- 
ing point; L is the latent heat of melting; AS m is the en- 
tropy variation at the melting point; dP/dT is the resulting 
Clausius-Clapeyron ratio at the melting point, a is the linear 
thermal expansion coefficient; RMSD is the root mean square 
displacement of atoms in the bulk solid at the melting point; 
S is the RMSD over the Na-Cl distance, for the Lindemann 
melting criterion. 
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FIG. 9: (a) Canonical evolution of the internal energy at 
liquid-solid coexistence. The decrease of internal energy at T 
= 1050 K indicates that the SL interface moves toward recrys- 
tallization. Conversely, the increase of internal energy at T — 
1070 K indicates that the liquid phase grows at the expenses of 
the solid phase. The internal energy for the BMHFT melting 
point of bulk NaCl lies in between, (b) Final microcanonical 
run showing the exact location of the melting temperature. 



The latent heat (enthalpy of melting) is estimated as 
the change of the internal energy at T m and found to 
be 0.29 eV/molecule. The calculated entropy jump at 
melting is thus AS m = L/T m = 6.32 ± 0.1 ks /molecule. 
A slope of the melting line (Clausius-Clapeyron relation) 
of 0.0311 kbar/K is then obtained. All 



dp 

dT 



these calculated results are in quite good agreement with 
experimental data. Tab. [n] summarizes the calculated 
values and their comparison with experiment. 



IV. CRYSTALLINE NACL(IOO): A 
NON-MELTING ANHARMONIC SURFACE 

Satisfied with the above description obtained for bulk 
NaCl, we moved on to study the NaCl surfaces, by sim- 
ulating slabs as described earlier. In slab simulations, 
we found the defect free NaCl(100) to warm up unevent- 
fully, and to remain solid and totally dry up to T m . The 
root mean square displacements (RMSD) of the first layer 
Na + and CD ions at T m were extracted and are shown 
in Fig. Elland Fig. El 
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FIG. 10: Solid NaCl instantaneous displacements of Na + ions 
in (a) bulk and (b) first surface layer at T m . 



Unsurprisingly, the surface ions vibrate more than bulk 
ions. The ratio of surface/bulk RMSD at T m is here 1.5. 
For comparison, in the Pietronero-Tosatti model [3^ the 
critical surface melting value of this ratio is ~ 1.6. More- 
over, vibrations are somewhat more extended in direction 
perpendicular to the surface than in the in-plane direc- 
tion. 

We subsequently verified that even above T m the 
NaCl(100) surface remained crystalline in a metastable 
state for at least 200 ps. In this overheated surface 
regime, solid NaCl(100) was found to possess a thick nu- 
cleation barrier against melting up to about T* = 1115 K 
~ T m + 50 K, in the following sense. When a thin surface 
film consisting of I atomic layers was artificially melted, 
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FIG. 11: Solid NaCl instantaneous displacements of Cl ions 
in (a) bulk and (b) first surface layer at T m . 
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and then let evolve to equilibrium at a grid of tempera- 
tures T > T m , the melted film was seen to spontaneously 
recrystallize for I < l cr it(T), with l cr it(T) > 1 for T < 
T* . (See Fig ll2|) . That indicates that in this tempera- 
ture range the overheated solid slab is locally stable, with 
a free energy barrier against overall melting Q . 

Above T*, overheating of crystalline NaCl(lOO) per- 
sisted until a higher "surface spinodal temperature" T ss ~ 
1215 K ~ T m + 150 K, now however with only a thin nu- 
clcation barrier £ cr jt(T) of a monolayer. 




1 — 1 — 1 — 1 — ■— ^ — 1 

1050 T M T 1150 T ss 1250 
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time. 

That explains why in practice we never even observed 
in simulation a spontaneous evaporation event off the 
solid surface. This however does not invalidate the sig- 
nificance of the simulation results. At any given instant 
of a realistically long time evolution of NaCl(100), there 
will be of course molecular evaporation and step flow, but 
that will still leave defect free terraces much larger on 
average than those we simulated. These terraces will be 
fully crystalline even well above T m , displaying precisely 
the microscopic non-melting behavior described above. 

These considerations however suggest that surface non- 
melting in alkali halides could not easily be pursued with 
static or nearly static experimental probes, because subli- 
mation and condensation will influence and possibly spoil 
the outcome. Perhaps the fast laser tools already em- 
ployed for metals |4(| could be brought to bear on this 
case too. The high vapor pressure also suggests tech- 
niques that do not rely on ultra-high vacuum [Zlj. We 
are currently considering in addition hard tip sliding fric- 
tion as a tool, with results to be described in forthcoming 
paper. 



V. SOLID SURFACE FREE ENERGY: 
THERMODYNAMIC INTEGRATION 



FIG. 12: Critical liquid film thickness vs. temperature for 
metastable solid NaCl(100) above T m . This corresponds to 
the effective thickness of the free energy barrier for the nucle- 
ation of the NaCl melt at the solid surface. 

Only at T ss , as high as 150 K above T m , does solid 
NaCl(100) become locally unstable and melt sponta- 
neously. We find the pronounced metastability of this 
solid surface to persist at least up to T* , or ~ 50 K above 
T m even in presence of common surface defects, such as 
molecular vacancies or steps. 

We have thus characterized NaCl(100) as a clear case 
of surface non-melting. This is a prediction that deserves 
to be tested in experiment. For a short enough time, it 
should be possible to raise the temperature of NaCl(100) 
by at least 50 K above T m , without any liquid sponta- 
neously nucleating at the surface. Of course, the solid 
surface will sublime very strongly at these temperatures. 
The experimental equilibrium vapor pressure of NaCl at 
T m is 0.345 mmHg, a large value. As a consequence, sur- 
face vacancies will form and surface steps will flow due 
to evaporation (and to re-condensation when working at 
solid- vapor equilibrium). We note that the time scale for 
this kind of surface evolution is a much slower one than 
that addressed here. The experimental rate of evapo- 
ration found from empirical expression [39^ valid for a 
variety of materials is about 4-10~ 5 gr/cm 2 s. For our 
typical 10x10 sized surface (200 surface NaCl units) of 
area around 3500 A 2 , and using an experimental vapor 
pressure of 0.345 mmHg the evaporation rate is ~ 1.7 
TO 5 s _1 . This indicates a typical evaporation time many 
orders of magnitude larger than our typical simulation 



In order to prepare for our thermodynamical assess- 
ment of the wetting of NaCl(100) by liquid NaCl, we 
eventually need to know according to eq. the solid- 
vapor interface free energy 7sv at the melting point. 

We calculated 7sv through thermodynamic integra- 
tion jlj using the following relation: 

where F is the free energy and E the internal energy. 
We simulated for this purpose a 2880 molecule bulk sys- 
tem and independently a slab comprising 1440 molecules 
and an equivalent volume of vacuum. We kept in this 
manner the all-important Ewald sum convergence un- 
changed, (we judge the implied extra error due to size 
effects to be negligible by comparison). By integrating 
the internal energy over (1/T) up to temperature T we 
separately obtained the bulk and the slab free energies 
per molecule up to T s and to T ss respectively. We did not 
explicitly include quantum freezing effects at low temper- 
atures, because they represent a small correction in com- 
parison with large thermal effects at the melting point. 
Nevertheless since quantum freezing takes place at tem- 
peratures T < (1/4)0d, (0d is the Debye temperature), 
we started integration at = (1/4)#d which is ~ 50 K 
for NaCl, therefore using the T = 50 K state as a refer- 
ence. 
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internal energies as a function of temperature. The slab con- 
sisted of 1440 NaCl molecules. 
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FIG. 14: The solid surface free energy from thermodynamic 
integration (line) and from the effective harmonic approxima- 
tion (dots). 



The bulk and slab internal energies as a function of tem- 
perature, are shown in Fig. 1131 After integration, the 
difference between slab and bulk free energies per unit 
surface area (accounting the presence of two surfaces) 
equals the surface energy. 

The NaCl(100) surface free energy calculated in this 
manner is displayed in Fig. El Starting at low tempera- 
ture, we note from the start a low surface energy, reflect- 
ing the excellent charge order of this surface. Upon in- 
creasing temperature, there is an important thermal drop 
of the surface free energy, especially fast for T > 600 K, 
indicating considerable additional anharmonicity of the 
solid surface relative to bulk solid NaCl. A good question 
is what part of that anharmonicity can be ascribed to 
temperature-dependent effectively harmonic vibrations, 
and what cannot. 

To judge on that, we computed separately the bulk 
and slab vibrational spectra as a function of tempera- 
ture. Using simulation trajectories in a bulk and in a slab 
comprising exactly the same number of molecules, we ex- 
tracted the Fourier transformed velocity autocorrelation 
functions of ions VNa(w) and Vci(w). The vibrational 
spectra are obtained as: 

%) = p|%M + ^%H- (13) 

k B l k B i 

By treating both the bulk and the slab spectra as a col- 
lection of harmonic oscillators, we obtained an effective 
surface vibrational free energy by subtracting their re- 
spective harmonic free energies off one another, and di- 
viding the outcome by two, for two surfaces. 

For increasing temperatures, the surface component of 
the slab spectra displays a stronger anharmonic softening 
than the bulk. This as anticipated gives rise to an "effec- 
tive harmonic" drop of surface free energy, as shown by 
dots in Fig. El We conclude that whereas about half the 
total anharmonic free energy decrease can be ascribed to 



the effective surface vibrational free energy and in partic- 
ular to the surface frequency softening with temperature, 
the remaining half cannot be accounted for in this way, 
representing "hard" anharmonicity. 

In conclusion the surface free energy of NaCl(100) ap- 
proximately halves its value from ~ 206 mJ/m 2 at 50 K 
to ~ 100 mJ/m 2 at T m . Such a large decrease of the 
already unusually small low-T surface energy results in 
an exceptionally stable solid surface. While the physical 
reason for a low surface energy at low-T is clearly the per- 
fect charge ordering, that for its large thermal decrease 
is the ability of rocksalt and of its surface to sustain ex- 
ceptionally large anharmonic vibrations without loss of 
mechanical stability. 



VI. THE LIQUID NACL SURFACE 

We proceeded next to study the liquid NaCl surface, 
or more correctly the liquid-vapor interface. There are 
earlier simulations of this liquid surface in the BMHFT 
model, notably by Heyes [lq. and there is theory work 
in the Restrictive Primitive Model 0j El • The more so- 
phisticated recent studies of Aguado et al. 01 of the liq- 
uid KI surface emphasize the role of polarization forces, 
and also summarize earlier work. 

Notwithstanding that, our scope here is to pursue a 
homogeneous comprehensive study of all surface proper- 
ties in a single simple model, and we therefore carried out 
a fresh study of NaCl in the BMHFT model, now with a 
larger size scale than that of Heyes [16| . 

Starting with the same 6x6x12 solid slab used in the 
previous section for the study of NaCl(100), temperature 
was raised above the surface spinodal temperature caus- 
ing the slab to melt. Subsequently, the liquid slab was 
gradually cooled down to T m , and equilibrated for about 
50 ps, after which correlations were examined. 
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Here too, evaporation of molecules and molecular 
dimers was very seldom observed, as an extremely rare 
event in our liquid surface simulations. The occasionally 
evaporated molecule traveled in vacuum to recondcnsc 
within a relatively short time, but had no chance to get 
otherwise equilibrated during the flight. The occasional 
molecular evaporation or condensation events were very 
quickly "forgotten" in the chaotic liquid surface dynam- 
ics. In this regime they therefore do not appear to affect 
at all the overall liquid surface behavior. This allows 
us to neglect our lack of an evaporation/condensation 
statistics, in that its inclusion would not alter the liquid 
surface properties to be extracted by the simulation. 

Because of the strong charge correlations in the bulk 
liquid, one might naively but not unreasonably have ex- 
pected this liquid surface to be structured, maybe lay- 
ered as in the metals Q, perhaps displaying a surface 
dipole ^ij . The actual liquid local surface density profile 

p H (z) obtained for both ionic species is shown in Fig. 1151 

All profiles are remarkably coincident and smooth^thus 
as one could also see from earlier MD studies to- 
tally devoid of layering. Moreover the Na and CI profiles 
are superposable with very great accuracy, thus the liquid 
surface displays no static average dipole either, whereas 
the local time and space dependent dipole fluctuations 
are large. 
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FIG. 15: Na, CI and total density profiles of the simulated 
liquid surface of NaCl at T = T m . Layering oscillations and 
surface dipoles are absent. Capillary fluctuations are very 
modest for our small cell size, and the large width of the 
liquid-vapor interface is a genuine local property of the sur- 
face. 

We worried that the large apparent surface spatial 
width should really represent the local liquid surface 
structure, and not, e.g., simply reflects long wavelength 
capillary fluctuations. For that purpose we carried out 
additional liquid slab simulations with alternatively a 
much smaller cell of size 4x4x4, or a much larger one 
of size 10x10x15. We found that the resulting surface 



density profiles remain essentially the same. We con- 
clude that the additional capillary broadening of the sur- 
face profile will only show up (logarithmically) for much 
larger sizes, and that the observed surface diffuseness is 
indeed intrinsic. 

The nature of diffuseness of the NaCl liquid surface is 
clarified by the simulation snapshot of Fig. 116b . show- 
ing very pronounced local thermal fluctuations in the in- 
stantaneous surface profile. This picture, suggestive of 
a low surface tension, high entropy surface, is in appar- 
ent contradiction with the massive non-wetting of solid 
NaCl(100) by its own melt, the latter implying a rela- 
tively high liquid surface tension. 



VAPOR 




FIG. 16: a) Simulation snapshot of the NaCl liquid surface 
at Tm. Note the large thermal fluctuations, with some nearly 
molecular configurations highlighted in the outermost region; 
b) Coordination number N(Z) and density profile, confirming 
a very smooth crossover from liquid (iV=4.6) to molecular 
vapor (iV=1.3 dotted line). 



In order to clarify the situation, we undertook a direct 
calculation of the surface free energy 7lVj equal to the 
surface tension, of liquid NaCl. The calculation was done 
by evaluating surface stress of the slab (which has two 
equivalent surfaces) via the Kirkwood-Buff formula: 
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(14) 



where: (a, /3) — (+, — ), ia and j/3 denote ions at site i or 
j, Z is the distance normal to the interface, L x , L y are the 
(x, y) dimensions of the supercell and <7ii = \{a xx + a yy ) 
and txj_ = a zz are the tangential and normal compo- 
nents of the stress tensor respectively. Here ( ) denotes 
a canonical average and Yin 1S over a ^ P arrs 01 P ar_ 



tides. Moreover r, 



Vij' 



is the interatomic 



distance, f a p(fij) is the force between atoms i and j, 
9a/3( r ij] Z) are the Na-Cl, Na-Na, Cl-Cl pair correlation 
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function measured in a slice centered at Z , p a (Z) the av- 
erage density of ion a near Z and finally A is a parameter 
here equal to one, but inserted for later use. 

The calculated liquid surface tension 7lv is shown as 
a function of temperature in Fig. 1201 The value at T m is 
104 ± 8 mJ/m 2 , in fairly good agreement with the exper- 
imental surface tension of 116 mJ/m 2 - and by chance 
essentially identical to that of the solid. A very large an- 
harmonicity was shown earlier to explain the relatively 
low surface free energy of solid NaCl(lOO). The physical 
reasons that make the liquid surface tension so relatively 
high will be addressed below. 

VII. LIQUID SURFACE ENTROPY DEFICIT: 
INCIPIENT MOLECULAR ORDER IN THE 
LIQUID SURFACE 

We wish to understand the reasons for the relatively 
high surface free energy of liquid NaCl. A clue is pro- 
vided by a comparison of solid and liquid surface excess 
entropies: 

S'surf = -d-y/dT. (15) 

Generally, one would expect that the much looser struc- 
ture and greater freedom of ionic motion at a liquid sur- 
face should yield a larger liquid surface entropy than that 
of the solid surface. Strikingly, in NaCl the smaller cal- 
culated temperature dependence of surface free energies 
shows just the reverse. We find a factor 2.6 lower surface 
entropy Slv compared with Ssv of the solid surface. Let 
us focus on this inverted result, which indicates quali- 
tatively speaking a liquid surface entropy deficit (SED). 
An entropy deficit is suggestive of some form of under- 
lying surface short range order. The order is clearly not 
layering: so what is it instead? 

The answer we found is that charge order, already very 
important in bulk, plays a newer and enhanced role at the 
molecular liquid surface. If surface thermal fluctuations 
are indeed very large, they are also revealingly correlated. 
For a Na + ion that instantaneously moves e.g., out of 
the surface, there is at least one accompanying CI - , also 
moving out; and vice versa. So while the large surface 
fluctuations smear the average liquid vapor density pro- 
file, bridging gently between the liquid and essentially 
zero in the vapor, (Fig. I15|) the two-body correlations, 
described e.g. by the the Na-Cl pair correlation function 

<7_| (r), and by its integral, the ion coordination number 

N, drop from values typical of the bulk liquid at T m to 
the nonzero value of the molecular vapor. For a quan- 
titative characterization, we calculated a locally defined 
charge coordination number: 

N ± (Z) = -L f +bZ fdZ'p(Z') [ d 3 r g + -(r;Z')\ 

... ^ 

where r m — 4.0A corresponds to the first local minimum 
of (r), and Sz is a small interval. Starting with Z 



inside the liquid slab, where the environment is bulk- 
like, we recovered N± L — 4.6 A at T m , as in the bulk 
liquid. Moving Z across the liquid-vapor interface we 
found N±(Z) to drop continuously from 4.6 downward 

(Fig. Hi)). 

In NaCl we know (even if simulation statistics is non- 
existing in the vapor) that N±{Z) for large Z is bounded 
below not by zero but by N± v , the average value for the 
NaCl vapor at T m . Experimentally the NaCl vapor con- 
sists for 69% of molecules (N=l), 31% of dimers, (N=2) 
and a trace of trimers |46j . The corresponding vapor av- 
erage is N± v ~ 1.3. Here emerges the crucial difference 
between the molecular NaCl vapor and e.g., the atomic 
LJ vapor, where N± v ~ 0. The larger the coordination 
number of ions in the surface region, the less their config- 
urational entropy, the higher the liquid surface tension. 
Hence incipient molecular order |45| could provide the 
reason for the SED found for liquid NaCl. 

For a test of this idea, we repeated the same Kirkwood- 
Buff calculation of the surface tension done previously, 
now however by slightly and artificially altering in 

Eq. (|14|) the value of correlations g^ at the surface. 

Specifically, we artificially reduced to zero the weight 
A attributed to forces acting among Na + and Cl~ ions 
for that (extremely small) fraction of outermost surface 
atoms whose coordination number N± < N± v ~ 1.3, 
the mean vapor value. The contribution of these config- 
urations to the pressure should provide a good measure 
of the the influence of incipient molecular charge order- 
ing to the surface free energy, in particular to the sur- 
face entropy. Through this highly artificial but in our 
view illuminating procedure, the surface internal energy 
(a mechanical variable) remains untouched, and thus only 
surface entropy is affected. 

We first identified the surface Na and CI atoms in the 
simulation by means of a simple algorithm. All ions are 
binned according to increasing Z and are represented by 
a sphere of finite radius (1.1 A in our case). An ion is 
considered a surface ion when the projection on the (x, y) 
plane of its representative sphere is non overlapping with 
that of any other ion at larger z. For each so identified 
surface ion i, we extracted from the simulation the in- 
stantaneous electrostatic potential value Vi, a quantity 
related to the coordination number, but more convenient 
to calculate and to handle. As shown in Fig.Elthe over- 
all electrostatic potential distribution of, say, Na ions in 
the liquid slab is shifted toward the electrostatic poten- 
tial value typical of the NaCl molecule, and away from 
that of the bulk liquid. This shift is evidently caused by 
the lower coordination of ions on the two surfaces of the 
slab, for the slab interior is bulk-like. The surface Na + 
ion potentials for example are shifted by ~ 0.1 eV on 
average relative to their bulk counterpart. The electro- 
static potentials of Cl~ ions behave specularly, and are 
thus shifted by ~ —0.1 eV relative to their bulk counter- 
parts. 

As the next step we established the necessary connec- 
tion between average electrostatic potentials (easily ex- 
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FIG. 17: The electrostatic potential distribution at Na + ions 
in (a) bulk liquid NaCl; (b) solid bulk, and (c) slab liquid 
of thickness 70 A, at T m . Comparison the potential at with 
BMHFT (d) trimers (e) dimers and (f) monomers is also pro- 
vided. Precise values of the average <f) are (a) —8.57 eV, (b) 
-8.32 eV, (c) -8.21 (d) -7.82 eV, (e) -7.485 eV, (f) -6.8 eV. 
The electrostatic distributions of CI - ions are just specular, 
that is identical to those of Na + , with a plus sign. 
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FIG. 18: Electrostatic potential distribution of surface Na + 
ions, plotted versus their own instantaneous coordination 
number. Superposed is the molecular correlation between po- 
tential and coordination(black dots). 



tracted from simulations, and at least in principle easily 
measurable) and coordination numbers (hard to extract 
from simulations, and probably harder to measure). The 
potential should vary monotonically with coordination, 
e.g., linearly for a fixed interatomic distance. Since in 
reality the Na-Cl distance increases with coordination, 
the overall dependence is somewhat less than linear as 
shown in Fig. El A raw histogram of electrostatic po- 
tentials against coordination numbers for the surface ions 
in the simulated liquid slab at T m is shown in Fig. 1181 

Finally, we computed a modified liquid surface tension 
by cutting off in the Kirkwood-Buff average in Eq. i|14|) 
the Coulomb part of the Na-Cl force contribution of 
surface Na ions, through a parameter A of the form 
A = 6(Vq — Vi) where 9 is the step function and Vq is 
a cutoff potential value, selecting the type of correlations 
to be removed. For example, Vb = —6.8 eV cuts off 
up to monomer correlations, Vq = —7.485 eV cuts off 
up to dimer correlations, etc. In particular we choose 
Vq = —6.99 eV, the value that cancels correlations for 
Na ions with TV < 1.3, the vapor average. 

With this tool, we are now able to examine more 
quantitatively the surface tension contribution due to 
the incipient surface molecular correlations causing 
N±(Z) — > N± v ~ 1.3, by correspondingly choosing the 
cutoff potential Vq — —6.99 eV. This modification gener- 
ally affects an exceedingly small fraction of surface ions. 
In particular, the Na - ions affected are quite few, as 
highlighted in Fig. 116b. Nevertheless the partial removal 
operated of the surface tension contribution due to this 
molecular part of surface Coulomb correlations yields a 
very considerable overall surface-tension decrease, with a 
large drop from 7lv = 104 mJ/m 2 to 7£ v = 53 mJ/m 2 , 



as shown in Fig. El and Fig. [201 This we interpret as a 




cp(eV) 

FIG. 19: The decrease of fictitious liquid surface tension 7l V 
at T 1T1 with increasing surface correlation coordination cutoff. 

direct confirmation that incipient molecular surface cor- 
relations are indeed responsible for the liquid SED and 
for the resulting high surface tension. Remarkably, since 
now 7£ v + 7sl < 7sV: the surface tension drop following 
the hypothetical removal of short range surface molecu- 
lar correlations would actually suffice to drive a complete 
instead of partial, wetting of NaCl(100) at the melting 
point. 

The increased temperature slope {dj^/dTl confirms 
that the calculated surface tension drop is directly re- 
lated to the restoring of a larger surface entropy, with re- 
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FIG. 20: Liquid NaCl surface tension. Note the good agree- 
ment of the calculated 7lv with experiment 7lv P . 7lv : arti- 
ficial surface tension calculated by setting A = for outer 
surface atoms with coordination below 1.3 (highlighted in 
Fig. \Wk ) . Once surface molecular order is removed in this 
way, surface entropy rises, surface tension drops, as shown. 
Solid NaCl(100) would be completely wet by this artificial 
liquid. 



moval of some of the SED through cancellation of molec- 
ular surface correlations. We found in fact that the drop 
from 7lv to 7^ v at T m corresponds exactly to the in- 
creased temperature slope —dj^y/dT, that is to the sur- 
face entropy increase, therefore with no change of sur- 
face internal energy as expected. In the presence of the 
surface molecular short-range order which we have just 
described, one could expect that the response to an ex- 
ternal electric field should be strongly reduced due to the 
effective neutralization. Indeed this effect is present and 
very visible in Heyes' early simulations of a BMHFT liq- 
uid slab in a parallel electric field, which further confirms 
our interpretation |l6j . 



VIII. SOLID-LIQUID INTERFACE FREE 
ENERGY: PARTIAL WETTING OF A NACL 
DROPLET ON NACL(IOO) 

In the above Sections we calculated the SV and the 
LV surface free energies of NaCl. The intermediate solid- 
liquid (SL) interface free energy 7sl, the third ingredient 
required to assess triple point wetting as in Eq. (JIJ is still 
missing. We calculated 7sl through a simulation of the 
partial wetting of solid NaCl(lOO) by a droplet of melt, 
and by using Young's equation Eq. © that connects it 
to the partial wetting angle |47j . 

A 500 molecule NaCl cluster was initially melted to 
form a nanodroplet. The droplet and the solid slab were 
separately equilibrated at 1050 K, and then brought to 
contact (fig. During the first 100 ps after contact, 



the droplet settled onto the substrate, slightly spreading 
and gradually approaching a final shape (fig. 121b). In the 
next 130 ps, spreading came to a halt, and the settled 
liquid droplet survived in a metastable, long lived state 
(Fig. I21b1. At the end of the simulation, the droplet- 
substrate system were as depicted in Fig. I21H . 
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FIG. 21: Time evolution of the liquid NaCl drop on 
NaCl(100). Dark and light circles stand for the Na + and 
Cl~ ions respectively. 

Let us consider the thermodynamics of this situa- 
tion. Because we are below T m (even if slightly) the 
true final equilibrium state should consist of a flat solid 
NaCl(100) surface, i.e. the nanodroplet should have com- 
pletely spread and recrystallized. That however will take 
a very long time. While the nanodroplet still exists, it 
forms an external wetting angle #lv (Fig. 12111 . as well as 
an internal angle 9$^. These angles obey the mechanical 
equilibrium equations Eq. J2J 9J. As it turns out, the 
angle $sl is irrelevant here, because it depends very crit- 
ically on temperature. In particular close enough to T m 
we expect 8sl — The external angle #lv is instead 
fully significant, as it depends very little on temperature, 
and indeed approximates the macro scop ic wetting angle 
measured in the bubble experiment [lj, Hj] ■ 

To determine the external wetting angle 6*lv of the 
nanodroplet we analyzed 100 configurations in a 100 ps 
equilibrated run. The instantaneous atomic positions 
were plotted in cylindrical coordinates (r and z, where 
r is parallel to the surface), and from the profile of the 
drop, we determined the best approximation to a spher- 
ical segment, by determining the center position and the 
radius. The contact angle follows immediately by simple 
geometry from these two quantities. Our best estimate 
obtained slightly below T m is 6*lv = 50° ± 5°. This value 
is in good agreement with the experimental value at the 
melting point (48°) [13j. At the end of the simulation 
the internal solid-liquid interface was still relatively sharp 
and flat, consistent with our assumption 0gL — 0. 

The connection between 6*lv and A7oo = (7lv +7sl) — 



14 



7sv is given directly by Young's equation: 

cos^v-1-^. (17) 

7LV 

When we plug in our calculated value of 7sv — 103 
mJ/m 2 , 7 Lv ~ 104 mJ/m 2 and finally 6 = 50°± 5°, we 
obtain 7sl = 36± 6 mJ/m 2 [48( and A7oo = 37 mJ/m 2 . 

A solid-liquid interface free energy of more than one- 
third the liquid surface tension is unusually large. Here, 
it is clearly attributable to an unusually large difference 
of density, as well as of correlations, between the solid 
and the liquid at the melting point. Corresponding to 
that, the solid- liquid interface is spatially rather abrupt, 
as shown by Fig. Inland Fig. [5] 

A further connection between A70C and the sur- 
face spinodal temperature T ss was made by assuming 
a phenomenological SL-LV interface interaction of the 
form [lOj: 

A700 s pLa - i) • (18) 

With our value of A7QQ, and L = 4.813 • 10 9 erg/g, a — 
5.9 A we predict T ss = 1210 K, quite close to that seen 
in simulations. 



IX. DISCUSSION AND CONCLUSIONS 

We studied in this paper the physics of the solid 
NaCl(lOO) surface and its wetting relationship with liq- 
uid NaCl at and near the melting temperature of bulk 
NaCl. Molecular dynamics simulations performed with 
classical BMHFT potentials were first of all shown to 
yield quite an accurate description of high temperature 
solid and liquid bulk. The NaCl(lOO) surface was sub- 
sequently studied, and found to be a non-melting sur- 
face, one that should in principle be possible to overheat 
well above T m . The solid surface free energy was cal- 
culated by thermodynamic integration and seen to drop 
very considerably at high temperature due to large an- 
harmonicities. The liquid NaCl surface was also stud- 
ied, and found to be very diffuse, strongly fluctuating, 
and devoid of static structure such as layering, or surface 
dipoles. However, calculation of the liquid surface ten- 
sion still gave a relatively large value, similar to the solid 
surface free energy at the melting point. The high surface 
tension signifies an unusual liquid surface entropy deficit, 
here ascribed to short range molecular order. In addition, 



the solid liquid interface free energy was also found to be 
relatively large, about one third the liquid surface ten- 
sion, consistent with the large density difference between 
solid and liquid. Direct simulation of a NaCl droplet de- 
posited onto NaCl(lOO) demonstrated very realistically 
the incomplete wetting, implied by the clear satisfaction 
of Eq. At T m , we obtained the solid, liquid and 

solid- liquid surface free energies, 103±4, 104±8mJ/m 2 , 
36±6mJ/m 2 respectively, quantitatively explaining the 
surface nonmelting of NaCl(lOO) through Eq. ©. We 
note that nonwetting of solid KC1 by its own liquid was 
arqued earlier by J. Rose and S. Berry ^| based on the 
behavior of a (KC1)32 clusters. 

The present work is meant as a prototype study, ide- 
ally repeatable for other elemental and molecular sys- 
tems including perhaps water and the surface of hexago- 
nal ice |50j . 

In alkali halides, the extraordinarily poor wetting of 
the solid (100) surface by the melt is traced to the con- 
spiracy of three separate factors, all of which can be fi- 
nally related to the long range Coulomb interaction be- 
tween ions: (i) surface anharmonicity stabilizes the solid 
surface; (ii) molecular correlations destabilize the liquid 
surface; (iii) a large density jump makes the solid-liquid 
interface very costly. 

Several aspects uncovered in, or implied by, this study 
should be amenable to direct experimental verification. 
X-ray studies should confirm the diffuseness of the liq- 
uid surface, the abruptness of the solid-liquid interface, 
and the molecular short range order at the liquid sur- 
face. Overheating of solid NaCl(100) should be directly 
observable. So might the temporary settling of liquid 
mini-droplets. Moreover, the much larger solid surface 
entropy should cause the partial wetting angle of liquid 
on solid NaCl to increase, rather then decrease with tem- 
perature. A brief account of some of the present results 
has appeared in Ref . [5lJ . 
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